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Numerical  discretization  for  unsteady  flow  simulations  can  be  broken  down  into  spatial  and  temporal  parts 
which  interplay  in  complex  and  sometimes  unexpected  ways.  This  paper  attempts  to  address  how  the  prop¬ 
erties  of  the  spatial  discretization  help  drive  the  choice  of  temporal  discretization.  In  addition,  it  examines 
methods  for  higher  than  second-order  accurate  time  integration  using  L-stable  singly-diagonally-implicit  (ES- 
DIRK)  Runge-Kutta  methods.1,2  Von  Neumann  analysis  is  used  to  examine  the  theoretical  effects  of  different 
spatial/temporal  discretization  combinations.  The  predictive  nature  of  the  von  Neumann  analysis  is  then  val¬ 
idated  through  the  exploration  of  the  convection  of  acoustic  waves  in  one  dimension  and  an  isentropic  vortex 
in  three  dimensions.  Is  is  shown  that  the  computational  results  follow  the  expected  trends  taking  the  von  Neu¬ 
mann  analysis  of  the  schemes  into  account.  This  work  highlights  that,  for  unsteady  problems,  both  dissipation 
and  dispersion  errors  must  be  accounted  for  when  selecting  optimal  Runge-Kutta  time  integrators. 


I.  Introduction 

The  use  of  high-order  spatial  discretizations  is  becoming  common  for  complex,  unsteady  flow  simulations.  How¬ 
ever,  the  accuracy  of  time  integration  methods  has  not  kept  pace  with  advances  in  spatial  accuracy  and  second-order 
accurate  methods  in  time  remain  the  norm.  This  limit  in  temporal  accuracy  is  primarily  a  result  of  the  second  Dahlquist 
barrier  which  states  that  no  A-stable,  implicit  linear  multistep  method  of  order  of  accuracy  greater  than  two  exists.*  1 * 
As  such,  for  higher-order  time  integration,  multistage  methods,  such  as  Runge-Kutta  (RK)  methods  must  be  utilized. 
Since  multistage  methods  require  multiple  residual  calculations  for  a  given  time  step,  they  seem  immediately  to  be 
poor  choices.  However,  Wang  and  Mavriplis  have  shown,  for  the  same  fourth-order,  diagonally  implicit  Runge-Kutta 
method  investigated  in  this  work,  that  high  accuracy  requirements  favor  higher-order  time  integration  using  larger  time 
steps  versus  lower-order  time  integration  with  smaller  time  steps.4  In  the  present  work,  an  attempt  is  made  to  gener¬ 
alize  these  results  to  a  broader  class  of  high-order  temporal  and  spatial  schemes.  Specifically,  von  Neumann  analysis 
is  performed  to  categorize  the  dissipation  and  dispersion  properties  of  several  candidate  schemes  and  to  discern  their 
comparative  strengths  and  weaknesses.  The  results  of  the  von  Neumann  analysis  are  then  shown  to  extend  readily  to 
real  world  flows  in  the  form  of  a  convecting  isentropic  vortex. 

Explicit  Runge-Kutta  time  integrators  may  seem  to  offer  a  cheaper  route  to  high-order  temporal  accuracy,  and, 
indeed,  are  widely  used  for  that  purpose.  However,  explicit  methods  are  limited  by  stability  restrictions  that  make 
them  unsuitable  for  stiff  problems  such  as  flows  with  disparate  physical  timescales.  Examples  include  low  Mach 
number  flows,  which  are  particularly  challenging  when  the  low-Mach  regions  occur  simultaneously  with  transonic  and 
supersonic  regions.  In  addition,  timescale  stiffness  occurs  in  high  Reynolds  number  boundary  layers,  wherein  the  close 
cell  spacing  in  the  wall-normal  direction  introduces  severe  time-step  restrictions  that  render  the  evolution  of  large-scale 
features  such  as  vortices  quite  inefficient.  Preconditioning  methods  have  been  developed  to  address  the  low-Mach 
limit,5,6,7,8’9’ 10,11,12  aspect  ratio  convergence  degradation,5- 8  and  both  the  low  and  high  Reynolds  number  limits,8,9 
among  others.  Moreover,  for  preconditioned  unsteady  solutions  in  physical  time,  a  dual-time  stepping  paradigm  must 
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be  used13  where  each  physical  time-step  is  treated  as  an  instantaneously  steady  solution  in  pseudo  time.  In  this  way, 
preconditioning  can  be  used  at  the  pseudo-time  level  and  the  accuracy  of  the  physical-time  solution  is  retained.  In  fact, 
it  has  also  been  shown  that  the  preconditioning  formulation  enables  the  definition  of  the  artificial  dissipation  terms  to 
be  cast  in  such  a  way  that  accuracy  is  preserved  for  stiff  problems.7  Therefore,  for  generality  and  because  of  interest  in 
multi-speed  and  high  Reynolds  number,  wall-bounded  flow  regimes,  a  dual-time  framework  is  adopted  in  the  present 
work. 

In  the  present  work,  von  Neumann  analysis14  is  used  systematically  to  examine  the  stability  behavior  and  the  nu¬ 
merical  dissipation  and  dispersion  errors  of  general  combinations  of  high-order  spatial  and  temporal  discretizations. 
Different  Runge-Kutta  time  integrators  are  applied  to  central-difference  spatial  schemes  with  added  artificial  dissipa¬ 
tion  terms.  The  overall  objective  is  to  gain  better  understanding  of  the  accuracy  of  these  schemes  and  to  determine 
optimal  choices  of  spatial  and  temporal  discretizations  for  practical  Euler  and  Navier-Stokes  simulations  at  all  speeds 
and  Reynolds  numbers. 

The  following  section  describes  the  governing  equations  and  spatial  and  temporal  discretizations.  Then,  the  theory 
behind  von  Neumann  analysis  as  well  as  its  results  in  the  form  of  dissipation  and  dispersion  error  analyses  are  pre¬ 
sented.  Next,  results  of  the  convection  of  acoustic  waves  in  one  dimension  and  an  isentropic  vortex  in  three  dimensions 
are  cataloged.  Finally,  important  conclusions  as  well  as  the  path  forward  from  the  present  work  are  summarized. 


II.  Governing  Equations 


Consider  the  Navier-Stokes  equations  written  using  dual-time  stepping  with  both  physical-  and  pseudo-time  steps, 
as  follows: 


3Q  5Q  3F, 

dx  dt  dxi  dxi  + 


(1) 


where  the  first  term  is  the  pseudo-time  derivative,  the  second  term  is  the  physical-time  derivative,  the  third  term  is  the 
convective  term  using  the  inviscid  fluxes  F(,  the  fourth  term  is  the  derivative  of  the  viscous  fluxes  V,,  and  the  last  term 
H  is  a  source  term.  In  the  present  work,  H  =  0  for  all  of  the  solutions  and  analyses  presented:  it  is  included  above 
for  thoroughness.  The  vector  of  conserved  variables  is  given  by  Q  =  [p  p»,  pt'o]7  ,  while  the  inviscid  fluxes  are 
given  by  F,-  =  [p iq  p u,Uj  +  p8ij  iiipho]7  where  h o  =  eo  +  p-  To  aid  in  analysis,  the  Navier-Stokes  equations  are 
frequently  written  in  quasi-linear  form  as  follows: 


5Q  9Q  3Q 

dt  ~dxi 


(2) 


where  A  =  S-  =  MAM  1 ,  with  M  and  M  1  being  the  right  and  left  eigenvectors,  respectively,  and  A  being  the  matrix 
of  eigenvalues,  X  =  {m,  +  c, n,  —  c}.  The  Navier-Stokes  equations  can  also  be  rendered  in  a  residual  formulation  as 


follows: 


3Q 

dx 


^  +  R,(Q)=o 


(3) 


where  Rv  =  ^  ^  —  H  is  the  complete  spatial  residual  of  the  chosen  spatial  discretization. 


A.  Spatial  Discretization 


In  the  present  work,  the  convection  term  will  be  discretized  using  central  differences  with  second,  fourth,  and  sixth 
orders  of  accuracy,  as  follows: 


3Tv  T,+i-T  ;■_! 

dxi  u  2Axj 


(4) 


3Ty  _  — Tj+2  + 810+1-810-! +Tj-2 

dxj  IV  12A Xj 


(5) 


3Tj_  Tj+3  -  9Y j+2  +  45Yj+i  -  45Tj_i  +  91,-2  -  Y;_3 

dxj  VI  60A Xj 


(6) 


where  T  could  be  F,  or  Q  depending  on  the  form  of  the  equations. 

As  is  well  known,  if  central  differences  are  used  for  the  convection  term,  it  is  generally  necessary  to  add  a  dissi¬ 
pative  term  to  the  Navier-Stokes  equations  to  ensure  convergence.  In  the  present  work,  scalar  artificial  dissipation  is 
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utilized  for  this  purpose  as  follows: 


„  _  3F d\i 
S  dxi  ^  11  X 11  dx?  dxf 


(7) 


where  r|  is  an  even  number  corresponding  to  one  more  than  the  order  of  accuracy  inherent  in  the  artificial  dissipation 
term,  ||  /,  1 1  =  \uj\  +  c  is  the  maximum  inviscid  eigenvalue  of  the  pseudo-linear  equations  where  ut  is  the  velocity  in  the 
direction  in  which  the  derivative  is  being  taken,  and  is  a  scaling  factor  appropriate  for  the  chosen  order  of  accuracy: 


£  TV  =  — 


Axf 

l2~ 


EVI  = 


Ax? 

60~ 


Typically,  is  discretized  using  a  second-order  accurate  finite  difference  formula  of  the  T| -degree  derivative.  The 

OXj 

order  of  accuracy  of  the  artificial  dissipation  term  is  given  by  the  power  of  the  Axj  in  the  coefficients. 


B.  Preconditioning 

To  accelerate  the  convergence  of  the  system  for  a  variety  of  physical  and  computational  phenomena  that  can  cause 
the  system  to  become  stiff,  preconditioning  can  be  used.  Phenomena  that  increase  the  stiffness  of  the  system  include, 
but  are  not  limited  to,  low  Mach  number,  high  aspect-ratio  cells,  low  and  high  Reynolds  number,  and  high-frequency 
unsteadiness.  Preconditioning  techniques  to  address  all  of  these  phenomena  have  already  been  developed  in  the  lit¬ 
erature.6'75'8'9  In  addition,  the  source  term  H  can  also  cause  the  system  to  become  stiff.  Preconditioning  to  address 
source  term  stiffness  could  also  be  developed  along  similar  lines,  depending  on  the  nature  of  the  source  term. 


C.  Temporal  Discretizations 


Thus  far,  the  focus  has  been  on  the  spatial  discretization.  In  this  section,  the  discretization  of  the  temporal  derivatives 
is  considered.  Both  physical-  and  pseudo-time  derivatives  are  discretized  with  Runge-Kutta  methods.  Consider  the 
flow  equations  written  in  residual  form  as  given  above,  but  repeated  for  convenience: 


dQ  dQ 

^  +  ^  =  -  R,(Q). 


First,  the  physical  time  derivative  will  be  discretized  as  follows: 


dQ  Q  "  -  Q" 
3x  +  At 


£amj  RS(Q7) 

j= 1 


(8) 


where  n  is  the  physical-time  step  index,  m  is  the  physical-time  Runge-Kutta  stage  index,  and  ,v  is  the  number  of  RK 
stages.  The  residual  weights  a,„j  are  drawn  from  the  Butcher  tableau  for  the  chosen  physical-time  method.  Each  stage 
occurs  at  a  discrete  fraction  of  the  time  step,  which  is  given  by  the  left  column  of  Cj  in  the  Butcher  tableau.  If  the 
Runge-Kutta  method  is  fully-implicit,  the  summation  in  equation  (8)  must  truly  range  over  all  stages  and  the  unknown 
stage  values  must  all  be  solved  simultaneously.  However,  if  a  diagonally-implicit  Runge-Kutta  method  is  used  instead, 
this  sum  proceeds  only  through  the  current  stage  because  the  Butcher  tableau  for  such  a  method  is  lower  triangular. 
Additionally,  since  current-stage  values  depend  only  on  past-stage  values,  only  the  current  stage  value  is  unknown 
during  each  stage  evaluation.  Thus,  subsequent  stages  can  be  solved  sequentially.  For  this  reason,  only  diagonally- 
implicit  Runge-Kutta  methods  are  considered  in  the  present  work.  The  updated  solution  at  the  end  of  the  physical-time 
step  is  calculated  from  the  already-solved  stage  values  as  follows: 

Q"  1 1  =  Q“  At  £  bjRs  (Q<)  (9) 

j=  i 


where  the  bj  are  taken  from  the  bottom  row  of  the  Butcher  tableau.  All  the  schemes  considered  herein  are  ’’stiffly 
accurate”  meaning  the  row  of  bj  is  identical  to  the  last  row  of  amj .  This  property  means  that  the  result  of  the  last 
stage  is  also  the  result  at  the  end  of  the  time  step  and  is  considered  to  be  an  essential  property  in  the  solution  of  stiff 
equations.11  The  Butcher  tableaux  for  all  of  the  Runge-Kutta  methods  mentioned  herein  are  given  in  Appendix  A. 

The  pseudo-time  derivative  can  be  discretized  either  implicitly  with  BDF1  or  explicitly  with  forward-Euler  or  an 
explicit  Runge-Kutta  method;  in  the  present  work,  Jameson’s  explicit  fourth-order  Runge-Kutta  is  used  for  pseudo¬ 
time  stepping.  Choice  of  the  pseudo-time  discretization  has  little  impact  on  solution  accuracy  and  primarily  impacts 
solution  efficiency.  As  such  the  precise  choice  of  pseudo-time  method  will  not  be  discussed  further. 
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III.  Von  Neumann  Analysis 


One  of  the  main  goals  of  this  analysis  is  to  examine  the  characteristics  of  different  time  integrators  when  applied 
to  specified  spatial  discretizations.  For  linear  equations,  von  Neumann  analysis  can  not  only  give  indications  of  the 
stability  of  a  given  combination  of  spatial  and  temporal  schemes,  it  can  also  indicate  how  the  numerical  solution  will 
differ  from  the  exact  solution.  The  output  of  von  Neumann  analysis  includes  both  of  the  following:  which  wave 
numbers  are  damped  and  by  how  much  per  time  step  (dissipation  error)  and  how  the  wave  speeds  at  various  wave 
numbers  correspond  to  the  exact  wave  speeds  at  those  same  wave  numbers  (dispersion  error).  The  following  sections 
describe  both  how  von  Neumann  analysis  can  be  applied  to  the  implicit  Runge-Kutta  time  integrators  and  the  results 
for  different  combinations  of  the  spatial  and  temporal  discretizations  discussed  above.  It  should  be  noted  that  although 
von  Neumann  analysis  is  technically  only  valid  for  linear  equations,  it  has  been  assumed  that  the  insights  garnered  from 
it  are  at  least  approximately  applicable  to  the  non-linear  Euler  equations.  It  should  also  be  noted  that  the  method  used 
in  this  paper  to  perform  von  Neumann  analysis  is  specifically  tailored  for  equations  that  can  be  treated  numerically  as 
systems  of  ordinary  differential  equations. 


A.  Theory 

To  derive  von  Neumann  analysis,  begin  with  the  scalar  equation: 


^7 

dt  dx 


Replace  q  with  its  Fourier  representation  at  spatial  points  j  and  temporal  points  n: 


„n  v00  ^co t  Jkx 

-  Lfc=i  Qke  e 

n  A  „(>>/  Jk(x-Ax) 

qj-l=Lk=l(ike  e  t  > 

tf"+i=E?=1fteaV*(*+ ^ 

qn+x  =  £“=1  qke^,+^eikx 

etc. 


(10) 


(11) 


Even  though  an  infinite  summation  is  shown  above,  a  given  resolution  can  in  fact  only  resolve  a  finite  number  of  modes. 
Taking  fa  =  1  (because  it  turns  out  that  all  the  fa  will  cancel  for  linear  systems)  and  using  the  above  substitutions  in 
equation  (10)  with  the  spatial  part  discretized  using  central  differences,  the  following  spatial  eigenvalue  contributions 
are  obtained  for  the  flux  contribution: 


Z(Q)n  =  --^XsinB 

Z(Q)iv  =  -  A.  [8  sin0 -sin  20] 

Z(B)vi  =  — 3^A.[45sin0  — 9sin20  +  sin30] 

Z(Q)viii  =  ~  i4oa.t^  [224 sin 0  —  56  sin 20+  y  sin  30  —  sin40] 

where  0  =  kAx  is  the  wave  number  which  ranges  from  [— n,  Jt].  For  the  artificial  dissipation  contribution,  the  following 
spatial  eigenvalue  contributions  are  obtained  for  second-order  accurate  r|-derivatives: 


Z(e)n  =  GnY(-l)5"1(2cos0-2)T1 


(13) 


where,  as  above,  r|  is  an  even  number  corresponding  to  one  more  than  the  order  of  accuracy  inherent  in  the  artificial 
dissipation  term  and: 


C// 


1 

6  Ax 


Cv/  = 


30Ax 


tvui 


1 

140Av 


which  mirror  the  coefficients  of  the  convective  terms  above.  It  should  be  noted  that  the  convective  contribution  is 
purely  imaginary  and  the  diffusive  contribution  is  purely  real.  The  total  spatial  eigenvalues  are  the  the  sum  of  the 
convective  term  for  the  order  of  accuracy  chosen  and  the  diffusive  term  for  the  order  of  accuracy  chosen. 

Once  the  eigenvalues  for  the  spatial  discretization  are  found,  they  are  then  used  in  the  time  discretization.  As  noted 
previously,  all  of  the  time  discretizations  considered  herein  can  be  written  as  Runge-Kutta  schemes.  The  amplification 
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factor  for  the  combined  spatial-temporal  scheme  can  be  calculated  as  follows  for  a  general  RK  scheme  with  a  Butcher 
tableau  as  referenced  in  the  previous  section: 


qn+l  =qll  +  AtYJbiki  (14) 

i=  1 

where  the  kj  are  the  stage  values  found  as  follows  for  equation  (10),  above: 

ki  =  Z(Q)qn+c'At  =Z(Q)  (qn+'paijk}j  .  (15) 

where  c,  is  the  fraction  of  the  time  step  at  which  the  current  stage  is  occurring,  as  given  by  the  left  column  of  a  Butcher 
tableau,  discussed  above.  By  rearranging  this  equation  into  the  following  form: 


(1  —aaZ)ki  —  z(  aijkj\=q"Z  (16) 

\j=h¥J  J 

it  becomes  clear  that  it  can  be  rewritten  as  a  linear  system  as  follows: 


1  —  a\\Z 

—auZ 

d\sZ 

k\ 

Z 

— ao\Z 

1  —  cmZ  • 

‘  d2sZ 

< 

ki 

>=qn< 

Z 

~as\Z 

a.aZ 

•  1  ClssZ 

ks 

z 

This  equation  is  then  solved  for  kj,  which  are  then  combined  in  equation  (14)  as  follows  to  find  the  amplification  factor 
G(Z)  as  a  function  of  the  complex  spatial  eigenvalues  as  follows: 


G 


=  At 


(18) 


When  the  magnitude  of  G(Z)  is  found,  |G|  gives  the  amplification  of  the  scheme  for  a  given  spatial  eigenvalue  or 
wave  number.  When  |G|  >  1,  the  scheme  is  unstable  as  the  values  are  amplified.  When  |G|  <  1,  the  scheme  is  stable 
as  the  values  are  dissipated.  Ideally  |G|  <  1,  but  also  very  close  to  one  for  all  wave  numbers  except  the  highest  wave 
number.  When  the  angle  formed  by  the  real  and  imaginary  parts  of  G  (Z)  is  found  and  divided  by  the  corresponding 
wave  speed  — A.0  for  the  given  scheme,  the  phase  error  of  the  overall  scheme  is  found,  as  follows: 


tan  1 


[r^g)  J 


-XQ 


(19) 


The  phase  error  §e  indicates  what  proportion  of  the  exact  wave  speed  the  wave,  corresponding  to  a  given  wave  number, 
is  traveling.  When  (])<,  >  1  the  wave  is  moving  faster  than  it  should  and  when  ©f,  <  1  the  wave  is  moving  slower  than  it 
should.  Ideally,  <|)f  =  1  at  all  wave  numbers. 


B.  Results 

The  results  of  the  von  Neumann  analyses  can  clarify  how  the  various  spatial  and  temporal  discretization  combinations 
behave.  For  instance,  when  progressively  higher-order  central  differences  are  used,  the  spatial  eigenvalues  can  be 
observed  to  stay  closer  to  their  exact  values  for  a  larger  range  of  wave  numbers.  This  trend  is  expressed  mainly  by 
shifting  the  location  at  which  the  extreme  values  occur  to  a  higher  wave  number  and  increasing  the  magnitude  of  that 
extreme,  as  can  be  seen  in  Figure  (1(a)).  In  other  words,  as  the  order  of  accuracy  increases,  the  eigenvalues  for  the 
lower  wave  numbers  become  more  and  more  linear.  Even  though  the  extreme  value  of  the  eigenvalues  increases  with 
increasing  order  of  accuracy,  the  condition  number  for  different  Mach  numbers  remains  the  same  no  matter  the  order 
of  accuracy  of  the  spatial  scheme  because  the  spatial  eigenvalues  scale  equally  for  all  wave  numbers. 

The  addition  of  artificial  dissipation  to  the  spatial  discretization  adds  a  negative  real  part  to  the  spatial  eigenvalues. 
As  the  order  of  accuracy  of  the  artificial  dissipation  term  increases,  the  magnitude  of  this  added  real  term  increases 
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Wave  Number  Wave  Number 


(a) 


(b) 


Figure  1 :  (a)  The  imaginary  part  of  the  spatial  eigenvalues  for  central  differences  of  second,  fourth,  sixth,  and  eighth 
order,  and  (b)  the  real  part  of  the  spatial  eigenvalues  for  artificial  dissipation  of  first,  third,  fifth,  and  seventh  order 


more  slowly  (indicating  better  accuracy),  as  can  be  seen  in  Figure  ( 1  (b)).  It  should  be  noted  that  since  scalar  dissipation 
is  used,  all  three  waves  experience  the  same  amount  of  dissipation  for  a  given  wave  number.  It  should  be  obvious  that 
this  can  cause  problems  at  low  Mach  numbers.  The  fastest  wave  has  a  ratio  of  dissipation  to  convection  of  one  at 
the  mid-wave  number  for  second-order  central  differences  and  first-order  artificial  dissipation,  whereas  the  slowest 
wave,  which  scales  like  Mach  number,  has  a  ratio  of  convection  to  dissipation  of  at  the  the  mid-wave  number.  It 
should  also  be  noted  that  the  CFL  number  acts  as  a  scaling  factor;  thus,  if  the  CFL  number  doubles,  both  the  real 
and  imaginary  parts  of  all  three  sets  of  spatial  eigenvalues  double.  For  this  reason,  it  is  sufficient  only  to  examine  the 
spatial  eigenvalue  trends  at  a  single  CFL  number. 

To  add  the  impact  of  the  temporal  discretization  to  the  von  Neumann  analysis,  the  complex  parametric  curves 
describing  the  spatial  eigenvalues  (with  the  wave  number  being  the  parameter)  are  overlaid  on  the  contour  map  of  the 
amplification  factor  in  the  complex  plane  for  the  chosen  temporal  scheme.  The  complex  valued  amplification  factors 
that  correspond  to  the  spatial  eigenvalue  curves  are  used  to  find  both  the  the  magnitude  of  the  growth  factor,  i.e.  one 
minus  the  numerical  dissipation  of  the  overall  scheme,  as  well  as  the  phase  error  of  the  overall  scheme,  as  a  function  of 
the  wave  number.  Figures  (2)  and  (3)  plot  the  magnitude  of  the  amplification  factor  versus  the  wave  number  for  BDF1 
and  ESDIRK4  temporal  discretizations.  In  all  cases,  only  the  most  quickly  propagating  wave  is  shown.  Subfigure  (a) 
charts  this  comparison  using  second-,  fourth-,  sixth-,  and  eighth-order  central  differences  with  no  artificial  dissipation, 
while  subfigure  (b)  plots  the  same  series  for  spatial  schemes  with  overall  accuracy  of  first  through  sixth  orders.  As 
can  be  seen,  somewhat  unexpectedly,  higher-order  central  differences  are  actually  more  dissipative  in  the  absence  of 
artificial  dissipation.  This  result  occurs  because  the  growth  factor  in  the  imaginary  direction  (y-axis)  decreases  fairly 
rapidly  for  BDF1 .  While  the  results  for  ESDIRK4  exhibit  this  same  trend,  it  is  important  to  note  the  scale  on  the  y-axis 
of  the  ESDIRK4  subplot  (a)  (i.e.  dissipation  is  negligible). 

Examining  subfigure  (b)  in  Figures  (2)  and  (3)  for  spatial  schemes  with  artificial  dissipation  and  effective  order 
of  accuracy  as  given  in  the  plot  legend,  it  is  clear  that  for  the  BDF1  scheme,  not  much  difference  exists  among  the 
spatial  discretizations  with  the  exception  of  the  first-order  accurate  discretization.  In  other  words,  the  BDF1  scheme 
itself  drives  the  overall  numerical  dissipation.  On  the  other  hand,  for  the  ESDIRK4  scheme,  it  is  clear  that  the  artificial 
dissipation  term  drives  the  overall  numerical  dissipation  of  the  scheme.  This  is  evidenced  by  the  fact  that  schemes 
with  overall  effective  orders  four  and  five,  for  instance,  which  share  a  common  artificial  dissipation  term,  are  visually 
coincident.  It  should  also  be  noted  that  almost  universally,  higher-order  spatial  discretizations  are  less  dissipative, 
especially  at  low  wave  numbers,  when  the  higher-order  ESDIRK4  time  integrator  is  used. 

Of  course,  the  dissipation  of  the  scheme  is  not  the  whole  story:  dispersion  error,  as  measured  by  the  phase  error, 
i.e.  how  fast  a  wave  travels  relative  to  how  fast  it  theoretically  should  be  traveling,  also  plays  a  key  role  as  a  source  of 
numerical  inaccuracy.  Figures  (4)  and  (5)  plot  the  phase  error  for  the  BDF1  and  ESDIRK4  time  discretizations  for  all 
the  same  spatial  discretizations  that  were  immediately  previously  discussed.  As  can  be  seen  for  the  backward-Euler 
time-discretization,  universally,  higher-order  central  differences  have  lower  phase  error  for  all  wave  numbers.  This 
trend  becomes  especially  apparent  when  the  overall  order-of-accuracy  discretizations  are  considered.  These  same 
observations  are  even  more  pronounced  when  the  ESDIRK4  time  discretization  is  considered  in  Figure  (5). 

Figure  (6)  plots  the  magnitude  of  the  growth  factor  (6(a))  and  the  phase  error  (6(b))  for  the  six  different  implicit 
time  integration  schemes  under  consideration  for  the  overall  fifth-order  accurate  spatial  discretization.  In  regards  to 
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Wave  Number  Wave  Number 


(a) 


(b) 


Figure  2:  The  magnitude  of  the  growth  factor  for  the  first-order  accurate  backward-Euler  time  discretization  at  a  CFL 
number  of  one  (a)  for  second-,  fourth-,  sixth-,  and  eighth-order  accurate  central  differences  with  no  artificial 
dissipation  and  (b)  for  the  combination  of  convection  and  artificial  dissipation  terms  that  produces  first  through  sixth 

spatial  orders  of  accuracy 


Wave  Number 


Wave  Number 


(a) 


(b) 


Figure  3:  The  magnitude  of  the  growth  factor  for  the  fourth-order  accurate  ESDIRK4  time  discretization  at  a  CFL 
number  of  one  (a)  for  second-,  fourth-,  sixth-,  and  eighth-order  accurate  central  differences  with  no  artificial 
dissipation  and  (b)  for  the  combination  of  convection  and  artificial  dissipation  terms  that  produces  first  through  sixth 

spatial  orders  of  accuracy 


Wave  Number 


Wave  Number 


(a) 


(b) 


Figure  4:  The  phase  error  for  the  first-order  accurate  backward-Euler  time  discretization  at  a  CFL  number  of  one  (a) 
for  second-,  fourth-,  sixth-,  and  eighth-order  accurate  central  differences  with  no  artificial  dissipation  and  (b)  for  the 
combination  of  convection  and  artificial  dissipation  terms  that  produces  first  through  sixth  spatial  orders  of  accuracy 
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Wave  Number  Wave  Number 


(a) 


(b) 


Figure  5:  The  phase  error  for  the  fourth-order  accurate  ESDIRK4  time  discretization  at  a  CFL  number  of  one  (a)  for 
second-,  fourth-,  sixth-,  and  eighth-order  accurate  central  differences  with  no  artificial  dissipation  and  (b)  for  the 
combination  of  convection  and  artificial  dissipation  terms  that  produces  first  through  sixth  spatial  orders  of  accuracy 


the  growth  factor,  all  of  the  ESDIRK  schemes  have  very  similar  growth  factor  curves,  especially  in  the  lower  half  of 
the  wave  numbers.  Crank-Nicolson  is  less  dissipative  than  any  of  the  three  aforementioned  schemes  at  all  except  the 
highest  wave  numbers,  likely  because  of  its  small  error  constant,  while  backward-Euler  is  more  dissipative  at  lower 
wave  numbers  and  less  dissipative  at  higher  wave  numbers,  which  is  not  a  desirable  arrangement.  When  it  comes 
to  phase  error,  generally,  more  accurate  temporal  discretizations  have  lower  phase  error,  as  can  be  seen  from  Figure 
(6(b)).  This  is  true  at  all  but  the  highest  wave  numbers.  The  exception  to  this  trend  is  BDF1,  which  has  the  highest 
phase  error  of  all  the  schemes  at  all  wave  numbers. 


Wave  Number  Wave  Number 


(a) 


(b) 


Figure  6:  The  (a)  magnitude  of  the  growth  factor  and  (b)  phase  error  for  the  time  discretizations  shown  in  the  legend 
at  a  CFL  number  of  one  with  overall  effective  fifth-order  spatial  accuracy 


Figure  (7)  plots  the  dissipation  (7(a))  and  dispersion  (7(b))  for  the  overall  fifth-order  accurate  spatial  discretization 
and  all  of  the  time  integrators  for  a  CFL  number  of  ten.  This  higher  CFL  number  is  plotted  to  assess  how  the 
spatial/temporal  schemes  respond  when  different  parts  of  a  grid  have  different  cell  sizes  (as  may  occur  in  the  presence 
of  boundary  layers)  and  to  ascertain  the  effect  of  using  larger  time  steps  if,  for  instance,  on  a  grid  that  would  produce 
the  desired  spatial  resolution,  the  spatial  error  is  still  the  dominant  form  of  error.  At  increased  CFL  numbers,  the 
schemes  are  much  more  differentiated,  with  Crank-Nicolson  clearly  being  the  least  dissipative  and  ESDIRK5  clearly 
being  the  least  dispersive.  In  light  of  the  fact  that  ESDIRK4  has  five  implicit  stages  whereas  ESDIRK5  has  seven 
implicit  stages,  the  potential  gains  of  ESDIRK5  over  ESDIRK4  do  not  necessarily  justify  the  increased  cost.  Indeed  in 
the  flow  results  section,  it  will  be  shown  that  for  the  problems  presented  in  this  work,  ESDIRK5  only  shows  worthwhile 
accuracy  gains  over  ESDIRK4  at  high  CFL  numbers  or  on  very  resolved  meshes. 

Overall,  the  von  Neumann  analysis  has  shown  that  when  a  high-order  spatial  discretization  is  used  a  high-order 
time  discretization  should  also  be  used.  Using  one  without  the  other  appears  to  be  an  inefficient  choice.  It  will  be  shown 
that  this  is  especially  true  on  more  refined  spatial  meshes.  This  analysis  also  shows  that  dispersion  and  dissipation  are 
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Figure  7:  The  (a)  magnitude  of  the  growth  factor  and  (b)  phase  error  for  the  time  discretizations  shown  in  the  legend 
at  a  CFL  number  of  ten  with  overall  effective  fifth-order  spatial  accuracy 


both  important  as  sources  of  error.  Indeed,  it  appears  that  dispersion  error,  in  general,  is  the  more  important  form  of 
error  since  higher-order  temporal  schemes  appear  mostly  to  correspond  to  schemes  with  lower  dispersion  error  but  not 
necessarily  lower  dissipation  error. 


IV.  Computational  Results 

Computational  results  are  presented  in  the  following  subsections.  First,  one-dimensional  solutions  of  the  Euler 
equations  are  presented.  Then,  results  of  the  three-dimensional  convection  of  a  two-dimensional  (cylindrical)  isen- 
tropic  vortex  are  considered. 


A.  One-dimensional  Results 


The  one-dimensional  Euler  equations  serve  as  a  simple  test  for  the  validation  of  the  dissipation  and  dispersion  char¬ 
acteristics  derived  in  the  von  Neumann  analysis.  The  periodic  governing  equations  utilize  conservative  variables  as 
expressed  in  equation  (1),  omitting  the  viscous  and  source  terms.  A  Mach  number  M„  =  0.5  is  chosen  and  the  uniform 
background  flow  is  debited  by  the  following  physical  parameters: 


p»„  =  8.7077  x  lO"1^,  ptioc  =  1.7458  x  102^,  =  400 K,  =  2.871  x  102j^,  y=  1.4  (20) 

While  the  Euler  equations  are  solved  in  divergence  form,  the  von  Neumann  analysis  assumes  linearity,  which  can  be 
accommodated  by  inducing  only  small  perturbations  with  magnitude  merely  one-percent  of  the  mean  bow.  The  bow 
is  perturbed  in  a  physically  consistent  manner  by  introducing  the  perturbations  using  the  characteristic  variables.  In 
this  way,  only  the  chosen  wave  numbers  0  =  kAx  can  be  excited.  This  process  is  given  by  the  following  equations: 

Qo  =  Q~  +  M8Qu.u±c  (21) 

5 Qu,u±c  =  8  ■  cos  (kx)  (22) 


where  8  =  0.01  or  one  percent  of  the  quiescent  bow  value,  as  stated  above.  The  subsequent  evolution  of  the  perturba¬ 
tion  wave  by  the  numerical  scheme  is  then  tracked  in  characteristic  variable  space,  so  that  it  can  be  directly  compared 
to  the  dissipation  and  dispersion  results  supplied  by  the  von  Neumann  analysis. 

Only  a  single  wavelength  L  of  the  perturbation  is  contained  in  the  domain.  This  wavelength  and  therefore  the 
corresponding  domain  length  is  set  to  one  meter.  Resolutions  I  of  ten  and  twenty  points  per  wave  are  used  so  it  can  be 
seen  that  there  are  necessarily  only  ten  and  twenty  points,  respectively,  in  the  domain.  These  resolutions  correspond 
to  a  wave  numbers  of  0  =  0.2ji  and  0  =  0.1 7t,  respectively.  Cumulative  dissipation  \G\cum  and  dispersion  error  tyecum 
after  a  number  of  time  steps  N  time  steps  can  be  predicted  from  the  von  Neumann  analysis,  as  given  by  the  following 
equations  (23): 


tyexum  — 


\G\cum 

Ay 

LXAcum 


\G\NU+C 

N-CFLu+c 

IL 


(i-«M 


(23) 
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where  Axcum  is  the  cumulative  absolute  distance  the  wave  is  out  of  phase.  Thus,  §e,cum  is  expressed  as  the  fraction  of 
the  wavelength  for  which  the  solution  is  out  of  phase.  The  following  computational  cases  focus  on  the  right  running 
acoustic  wave  ( qu+c )  and  are  compared  to  theoretical  results  from  von  Neumann  analysis. 

Table  (1)  and  Figure  (8)  compare  the  performance  of  Crank-Nicolson  and  the  three  ESD1RK  schemes  using  an 
overall  fifth-order  accurate  spatial  discretization  with  artificial  dissipation.  There  are  ten  points  per  wave,  a  CFL 
number  of  unity  is  used,  and  the  wave  is  advanced  one  period.  As  can  be  seen,  all  schemes  exhibit  low  dissipation, 
while  Crank-Nicolson  exhibits  noticeable  dispersion.  Examining  the  error  plot,  it  becomes  clear  that  the  higher  the 
order  of  the  scheme,  the  lower  error  it  exhibits,  with  very  little  difference  apparent  between  ESD1RK4  and  ESDIRK5. 
The  following  pattern  should  also  be  observed  and  noted:  when  the  error  is  primarily  dispersive,  the  minimum  error 
will  occur  at  a  position  near  the  extrema  of  the  waveform;  however,  when  the  error  is  primarily  dissipative,  the 
minimum  error  will  occur  at  a  position  near  the  inflection  points  of  the  waveform. 

Table  1 :  Von  Neumann  and  numerical  results  for  the  forward  moving  acoustic  wave  for  the  Euler  equations  with 
artificial  dissipation,  with  ten  points  per  wave,  at  a  CFL  number  of  one,  after  one  period  of  convection 


Dissipation  Error 

Dispersion  Error 

Scheme 

VNA 

Simulation 

VNA 

Simulation 

CN 

8.42  x  10“3 

2.77  x  1(T2 

3.14  x  1(T2 

3.08  x  10“2 

ESDIRK3 

4.25  x  10~2 

4.26  x  10~2 

2.47  x  1(T3 

2.39  x  10‘2 

ESDIRK4 

9.25  x  10“3 

9.26  x  10“3 

5.37  x  1(T4 

5.26  x  10~4 

ESDIRK5 

9.31  x  10“3 

9.31  x  10“3 

4.10  x  10~4 

4.00  x  10~4 

X-Coordinate 


(a) 


(b) 


Figure  8:  (a)  Amplitude  and  (b)  error  in  the  convection  of  the  forward  moving  acoustic  wave  for  the  Euler  equations 
including  artificial  dissipation,  with  ten  points  per  wave,  at  a  CFL  number  of  one,  after  one  period  of  convection 


Table  (2)  and  Figure  (9)  shows  the  same  data  as  previously  for  the  same  CFL  number  and  resolution,  but  after  ten 
periods  of  convection.  As  is  clearly  seen,  error  accumulates  as  solution  time  increases,  with  Crank-Nicolson  exhibiting 
much  more  dispersion  error  and  ESDIRK3  showing  much  more  dissipation  error.  ESDIRK4  and  ESD1RK5  even  show 
a  small  amount  of  dissipation,  but  the  differences  between  them  are  still  mostly  unnoticeable. 

Table  2:  Von  Neumann  and  numerical  results  for  the  forward  moving  acoustic  wave  for  the  Euler  equations  with 
artificial  dissipation,  with  ten  points  per  wave,  at  a  CFL  number  of  one,  after  ten  periods  of  convection 


Dissipation  Error 

Dispersion  Error 

Scheme 

VNA 

Simulation 

VNA 

Simulation 

CN 

8.11  x  10~2 

8.49  x  10~2 

3.14  x  10-1 

3.14  x  10"1 

ESDIRK3 

3.52  x  10'1 

3.60  x  lO'1 

2.47  x  10~2 

2.41  x  10‘2 

ESDIRK4 

8.88  x  10~2 

8.93  x  10”1 

5.37  x  10~3 

5.22  x  10‘3 

ESDIRK5 

8.93  x  10~2 

8.96  x  10~2 

4.10  x  10~3 

4.00  x  10‘3 
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(a)  (b) 


Figure  9:  (a)  Amplitude  and  (b)  error  in  the  convection  of  the  forward  moving  acoustic  wave  for  the  Euler  equations 
including  artificial  dissipation,  with  ten  points  per  wave,  at  a  CFL  number  of  one,  after  ten  periods  of  convection 


For  all  of  the  remaining  one-dimensional  cases,  resolution  has  been  increased  to  twenty  points  per  wave.  Table  (3) 
and  Figure  (10)  chart  a  case  with  the  same  CFL  number  of  one  and  ten  periods  of  progression,  as  previously,  but  at  this 
increased  spatial  resolution.  The  dispersive  properties  of  Crank-Nicolson  are  still  evident  at  this  finer  resolution,  but 
the  dissipation  of  ESDIRK3  appears  to  have  been  mostly  eliminated.  On  examination  of  the  error  subplot,  it  becomes 
clear  ESDIRK3  still  exhibits  much  higher  dissipation  than  either  ESDIRK4  or  ESDIRK5. 

Table  3:  Von  Neumann  and  numerical  results  for  the  forward  moving  acoustic  wave  for  the  Euler  equations  with 
artificial  dissipation,  with  twenty  points  per  wave,  at  a  CFL  number  of  one,  after  ten  periods  of  convection 


Dissipation  Error 

Dispersion  Error 

Scheme 

VNA 

Simulation 

VNA 

Simulation 

CN 

3.05  x  10~3 

1.00  x  10“2 

8.11  x  10~2 

8.11  x  10~2 

ESDIRK3 

5.02  x  10“2 

5.02  x  10'1 

1.51  x  10~3 

1.53  x  10“3 

ESDIRK4 

3.13  x  10~3 

3.13  x  10“3 

1.50  x  10~4 

1.58  x  10~4 

ESDIRK5 

3.14  x  10~3 

3.14  x  10“3 

6.78  x  10-5 

6.90  x  10“5 

Implicit  time  integration  methods,  like  the  ones  presented  herein,  are  unconditionally  stable  and  as  such  can  be 
run  at  higher  CFL  numbers  than  stability  considerations  would  allow  for  explicit  schemes.  Increased  CFL  numbers 
result  in  larger  time  steps,  which  increases  the  amount  of  temporal  error.  If  the  solution  is  already  in  a  regime  where 
spatial  error  dominates,  an  increased  CFL  number  will  not  greatly  affect  the  overall  error;  however,  if  temporal  error  is 
dominant,  the  overall  accuracy  of  the  solution  will  be  degraded.  Results  for  simulations  at  CFLu+c  =  10  and  which  are 
advanced  for  only  two  time  steps,  corresponding  to  one  solution  period,  are  shown  in  Table  (4)  and  Figure  (11).  In  this 
figure,  the  differences  in  the  accuracy  among  the  various  temporal  schemes  under  consideration  are  clear.  As  the  von 
Neumann  analysis  predicted.  Crank  Nicolson  preserves  amplitude  well  yet  is  too  dispersive,  with  its  solution  lagging 
about  36%  of  the  wavelength.  ESDIRK3  exhibits  high  dissipation  and  dispersion,  damping  50%  of  the  amplitude  and 
lagging  the  exact  solution  by  20%  of  the  wavelength.  ESDIRK4  and  ESDIRK5  are  still  very  well  behaved  with  the 
solutions  using  them  showing  99%  and  95%  of  the  original  amplitude  and  phase  errors  of  only  about  5%  and  1%  of 
the  wavelength,  respectively).  However,  the  superior  precision  of  ESDIRK5  is  evident  for  this  case;  supporting  the 
predictions  of  von  Neumann  analysis,  the  computational  results  show  that  the  fourth-order  ESDIRK  scheme  is  less 
dissipative  yet  more  dispersive  than  the  fifth-order  scheme. 

Finally,  Table  (5)  and  Figure  (12)  performance  of  the  temporal  schemes  under  consideration  for  a  very  long  solution 
time.  In  fact,  each  wave  is  advanced  for  twenty-thousand  time  steps  at  a  CFLu+c  =  1,  meaning  that  the  wave  convects 
for  one-thousand  wave  lengths.  The  popular  Crank-Nicolson  scheme  does  well  to  preserve  the  magnitude  of  the 
wave,  which  has  74%  of  its  original  amplitude  at  the  end  of  the  computation.  However,  the  Crank-Nicolson  solution 
suffers  tremendously  from  dispersion  error.  To  best  capture  the  wave  features,  the  aforementioned  figure  features 
only  two  wavelengths,  which  obscures  the  fact  that  the  Crank-Nicolson  solution  is  actually  about  eight  wavelengths 
behind  its  correct  position.  As  the  plot  shows,  ESDIRK3  is  obviously  extremely  dissipative  and  effectively  damps 
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(a)  (b) 


Figure  10:  (a)  Amplitude  and  (b)  error  in  the  convection  of  the  forward  moving  acoustic  wave  for  the  Euler  equations 
including  artificial  dissipation,  with  twenty  points  per  wave,  at  a  CFL  number  of  one,  after  ten  periods  of  convection 


Table  4:  Von  Neumann  and  numerical  results  for  the  forward  moving  acoustic  wave  for  the  Euler  equations  with 
artificial  dissipation,  with  twenty  points  per  wave,  at  a  CFL  number  of  ten,  after  one  period  of  convection 


Dissipation  Error 

Dispersion  Error 

Scheme 

VNA 

Simulation 

VNA 

Simulation 

Crank-Nicolson 

0.01  x  10-' 

0.02  x  10'1 

3.61  x  lO"1 

3.50  x  lO'1 

ESDIRK3 

4.90  x  10-' 

4.90  x  10'1 

1.92  x  10"1 

2.00  x  10'1 

ESDIRK4 

0.07  x  10-1 

0.07  x  10'1 

4.90  x  10~2 

5.00  x  10~2 

ESDIRK5 

0.51  x  10-' 

0.55  x  10”1 

1.39  x  10“2 

1.00  x  10~2 

(a)  (b) 


Figure  11:  (a)  Amplitude  and  (b)  error  in  the  convection  of  the  forward  moving  acoustic  wave  for  the  Euler  equations 
including  artificial  dissipation,  with  twenty  points  per  wave,  at  a  CFL  number  of  ten,  after  one  period  of  convection 
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the  perturbation  entirely.  ESD1RK4  and  ESDIRK5,  on  the  other  hand,  exhibit  excellent  dissipation  and  dispersion 
characteristics,  especially  given  the  long  solution  time.  For  both  schemes,  damping  is  kept  to  approximately  27%  as 
the  final  solution  maintains  about  73%  of  the  the  correct  amplitude.  Dispersion  is  «  3%  of  the  wavelength.  The  ability 
of  ESDIRK4  and  ESDIRK5  to  preserve  amplitude  and  phase  in  the  presence  of  artificial  dissipation  seems  highly 
attractive  for  time-accurate  solutions  with  high-order  spatial  discretizations.  Figuref  12(b))  shows  the  exact  solution 
along  with  the  absolute  error  of  ESDIRK4  and  ESDIRK5  schemes,  capturing  the  combined  effect  of  cumulative 
dissipation  and  dispersion  errors  and  highlighting  that  ESDIRK5  does,  in  fact,  have  slightly  lower  overall  error  than 
ESDIRK4. 

Table  5:  Von  Neumann  and  numerical  results  for  the  forward  moving  acoustic  wave  for  the  Euler  equations  with 
artificial  dissipation,  with  twenty  points  per  wave,  at  a  CFL  number  of  one,  after  one-thousand  periods  of  convection 


Dissipation  Error 

Dispersion  Error 

Scheme 

VNA 

Simulation 

VNA 

Simulation 

CN 

2.63  x  KT1 

2.65  x  10 1 

8.11  x  10° 

8.10  x  10° 

ESD1RK3 

9.94  x  KT1 

9.94  x  10'1 

1.51  x  lO^1 

1.00  x  10"1 

ESDIRK4 

2.69  x  KT1 

1.95  x  10'1 

1.50  x  10~2 

3.00  x  10‘2 

ESDIRK5 

2.70  x  10'1 

2.01  x  10'1 

6.78  x  10~3 

2.50  x  10‘2 

(a)  (b) 

Figure  12:  (a)  Amplitude  and  (b)  error  in  the  convection  of  the  forward  moving  acoustic  wave  for  the  Euler  equations 
including  artificial  dissipation,  with  twenty  points  per  wave,  at  a  CFL  number  of  one,  after  one-thousand  periods  of 

convection 

These  one-dimensional  Euler  computations  serve  as  a  validation  of  the  von  Neumann  analysis  results  and  have 
helped  inform  perceptions  of  the  strengths  and  weaknesses  of  the  temporal  schemes  considered  herein.  Having  estab¬ 
lished  expected  performance  near  the  linear  regime  of  the  Euler  equations,  the  next  section  presents  a  more  complex, 
three-dimensional  case. 

B.  Convection  of  an  Isentropic  Vortex 

The  isentropic  vortex  is  a  generally  localized  set  of  specific  perturbations,  which  occur  across  a  broad  range  of  wave 
numbers,  to  an  otherwise  uniform  flow.  The  perturbations  occur  in  the  x-  and  v-velocity  components  and  the  tempera¬ 
ture.  Because  these  perturbations  are  isentropic,  they  result  in  perturbations  in  the  pressure,  density,  and  energy  of  the 
flow  as  well.  The  velocity  and  temperature  perturbations  of  the  isentropic  vortex  are  as  follows: 


5m  =  -VRccToc^-  (y-yo)e^1  '  ) 

2n 

(24) 

Sv  =  y/RooToo-^-  (x  —  xo)e^(l~r2') 

2n 

(25) 

(26) 

16([>y k- 
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where  a  determines  the  strength  of  the  vortex  and  <f>  determines  the  gradient  of  the  vortex,  which  in  turn  determines 
vortex  size.  Higher  values  of  <|>  correspond  to  smaller  vortices  with  more  rapidly  changing  temperature  and  velocity 
from  the  edge  to  the  center.  The  vortex  center  is  at  (jco,yo)  with  the  distance  to  this  center  from  a  coordinate  (x,y) 

given  as  r  =  \J (x  —  xo)2  +  (y  —  yo)2.  If  these  perturbations  look  a  little  different  from  those  in  other  publications,  it 
is  because  they  are  in  dimensional  form.  Many  other  publications  assume  a  particular  non-dimensionalization,  but  do 
not  always  explicitly  state  that  fact.4 

For  all  of  the  vortex  test  cases  presented  herein  a  =  4  and  <f»  =  1 ;  as  a  result,  the  vortex  has  a  diameter  of  about  2.2 m 
(vortex  diameter  is  determined  as  the  diameter  of  the  area  where  the  density  perturbation  is  more  than  five-percent  of 
the  free  stream  value).  The  cases  are  run  at  a  free  stream  Mach  number  of  M„  =  0.5  in  the  x-direction.  The  free  stream 
density  is  set  to  be  poo  =  1.0^-.  It  was  desired  that  the  free  stream  speed  of  sound  be  a  round  number  close  to  its  typical 
value  for  air,  so  c„,  =  400y  was  chosen.  Given  that  R„,  =  287. 1 1  and  y  =  1 .4  for  air,  this  necessarily  means  that 
free  stream  temperature  and  pressure  have  the  following  values  7'x,  =  398. 06 A"  and  p„  =  1 14285.7/J«.  Thus,  the  free 
stream  conserved  variables  have  the  following  values: 


,  =  1 .0— 


m3  7 


pUoo  =  200.0- 

r. 


pVoo  =  0.0- 

r, 


pvi’oo  =  0.0- 

r, 


pe0.oo  =  305714.3- 


(27) 


Given  the  strength  of  the  vortex  and  the  free  stream  Mach  number,  this  test  case  will  clearly  be  affected  by  the  non¬ 
linearity  of  the  Euler  equations 

First,  the  formal  orders  of  accuracy  of  both  the  spatial  and  temporal  discretizations  examined  herein  are  establishes 
by  convecting  the  isentropic  vortex,  described  above,  for  a  short  amount  of  time.  For  spatial  accuracy,  the  vortex  is 
convected  for  one  time  step  whose  fixed  length  is  At  =  0.00001 0742 1875s  on  a  computational  mesh  that  is  22.0m  x 
22.0 m  in  the  x-  and  y-directions  and  whose  resolution  varies.  Because  the  flow  should  have  no  variation  in  the  z- 
direction,  four  points  were  used  in  that  dimension  for  all  grid  resolutions.  As  a  result,  the  length  of  the  grid  in  z  varies 
depending  on  the  x  and  y  grid  resolution  such  that  Ay  =  Ay  =  A z.  Thus,  for  a  grid  resolution  of  one-hundred  nodes 
in  the  x-  and  y-directions,  the  mesh  would  be  0.88 m  in  the  z-direction.  Boundaries  are  periodic  in  all  three  directions. 
The  given  time-step  size  corresponds  to  a  CFL  number  close  to  one  when  the  mesh  resolution  is  1600  points  in  each 
of  the  v-  and  y-directions.  Figure  (13)  plots  the  results  of  the  spatial-order-of-accuracy  tests.  As  can  be  seen,  the  errors 
of  all  the  spatial  schemes  considered  exhibit  the  expected  convergence  trends  as  resolution  increases.  Since  both  axes 
are  given  in  terms  of  powers  of  two  in  Figure  (13),  a  third-order  scheme,  for  example,  is  expected  to  decrease  three 
grid  lines  for  every  one  grid  line  decrease  in  dx. 


Figure  13:  Variation  of  the  Z.2-norm  of  solution  error  in  density  with  spatial  resolution  for  the  isentropic  vortex  after 
one  time  step,  demonstrating  formal  order-of-accuracy  convergence  for  the  given  spatial  schemes 


For  the  temporal  accuracy  study  the  vortex  is  convected  for  0.0064.?  on  a  computational  mesh  that  is  24.0 m  x 
24.0m  x  0.04m  with  a  fixed  resolution  of  2400  x  2400  x  4  nodes  in  each  of  the  three  directions,  yielding  Ay  =  Ay  = 
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A z  =  0.01.  Again,  boundaries  are  periodic  in  all  three  dimensions.  This  case  is  run  with  different  time  step  sizes 
ranging  from  as  large  as  0.0064s-  to  as  small  as  6.25  x  1 0  6.v  with  each  smaller  time  step  being  j  as  large  as  its 
immediately  bigger  time  step.  These  time  steps  correspond  to  CFL  numbers  ranging  from  896  to  0.875.  Obviously, 
this  means  that  for  smaller  time-step  sizes,  the  solution  must  be  run  for  a  greater  number  of  time  steps  to  achieve  the 
overall  0.0Q64.V  of  solution  time.  Figure  (14)  plots  the  density  error  convergence  of  Crank-Nicolson  as  well  as  third-, 
fourth-,  and  fifth-order  ESDIRK.  As  can  be  seen,  all  four  schemes  demonstrate  the  expected  order-of-accuracy  error 
convergence. 


Figure  14:  Variation  of  the  Z.2-norm  of  solution  error  in  density  with  time-step  size  for  the  isentropic  vortex  after 
0.0064  seconds,  demonstrating  formal  order-of-accuracy  convergence  for  the  given  temporal  schemes 

Next,  the  vortex  with  conditions  as  given  above,  is  convected  for  a  much  longer  time  on  the  same  mesh  that  was 
used  for  the  spatial-accuracy  study,  i.e.  22.0m  x  22.0m  in  the  x-  and  y-directions.  Since,  as  was  stated  earlier,  the 
vortex  diameter  is  about  2.2 m,  the  domain  is  about  ten  vortex  diameters  long  and  there  are  eleven  points  across  the 
vortex  on  a  mesh  with  one-hundred  nodes  in  x  and  y.  The  density  perturbation  of  the  vortex  along  its  x-direction 
center-line  after  two  and  four  domain  lengths  of  convection  using  the  stated  grid,  at  a  free  stream  CFL  number  of 
about  one  (1.09375  to  be  precise),  and  for  the  four  different  time  discretizations  is  shown  in  Figure  (15).  As  can  be 
seen,  all  three  ESDIRK  schemes  produce  the  same  results  while  Crank-Nicolson  shows  only  a  small  variation,  even 
after  forty  vortex  widths  of  convection.  This  results  seems  to  indicate  that  for  this  mesh  at  a  CFL  number  of  one, 
spatial  error  dominates  temporal  error  for  all  time  schemes  under  consideration. 

To  verify  the  spatial-error-dominance  hypothesis.  Figure  (16(a))  plots  solution  error  after  a  short  time  (0.02754) 
versus  grid  resolution  at  a  constant  CFL  number  of  about  one.  Because  the  CFL  number  is  held  constant,  the  temporal 
resolution  necessarily  increases  at  the  same  rate  as  the  spatial  resolution  from  right  to  left  along  the  x-axis  of  the  chart. 
As  can  be  seen,  all  three  ESDIRK  schemes  as  well  as  the  explicit  third-order  Runge-Kutta-Wray,  which  is  frequently 
used  as  the  time-integrator  in  high-order  Cartesian  codes  and  is  included  for  the  sake  of  comparison,  have  fifth-order 
error  convergence.  Crank-Nicolson  initially  shows  fifth-order  error  convergence,  but  its  accuracy  reduces  to  second- 
order  as  resolution  increases,  signaling  that  at  higher  resolution,  the  temporal  scheme  becomes  the  dominant  source 
of  error.  It  should  be  noted  that  the  one-hundred  node  grid  used  above  corresponds  to  the  second  point  from  the  right 
(as  noted),  meaning  it  is  clearly  in  the  spatially-dominant  regime  for  all  time  schemes.  Figure  (16(b))  shows  the  error 
of  the  different  time  schemes  relative  to  the  the  ESD1RK5  scheme  at  a  CFL  number  of  about  one  (to  allow  for  more 
straightforward  visualization  of  the  differences).  As  can  be  seen,  at  finer  resolutions,  ESDIRK3  has  about  half  the 
error  of  RK-Wray,  while  ESDIRK4  and  ESDIRK5  have  the  same  amount  of  error  for  the  entire  domain  of  the  plot, 
signaling  that  spatial  error  is  clearly  dominant  for  these  time  schemes  at  all  tested  grid  resolutions. 

Referring  again  to  Figure  (16),  it  can  be  seen  that  Crank-Nicolson  starts  to  have  a  different  slope  at  the  third  point 
from  the  left,  which  corresponds  to  forty-one  points  across  the  vortex,  or  four  times  the  resolution  shown  previously. 
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X-Coordinate 


X-Coordinate 


(a) 


(b) 


Figure  15:  Density  at  the  vortex  center  line  for  the  time  discretization  given  in  the  legend  with  eleven  points  across 
the  vortex  at  a  CFL  number  close  to  one  after  convecting  (a)  two  periods,  (b)  four  periods  . 
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Figure  16:  Comparison  of  the  density  error  for  different  time  discretizations  at  a  CFL  number  of  one  for  (a)  absolute 
error  and  (b)  error  scaled  by  the  ESD1RK5  error  at  a  CFL  number  close  to  one. 
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A  significant  difference  in  error,  however,  is  not  seen  until  the  next  point  to  the  left  with  eighty-one  points  across  the 
vortex.  Figure  (17)  shows  the  density  at  the  vortex  center  line  for  the  vortex  on  this  eight-times  resolution  grid  after 
four  periods  of  convection.  All  lines  appear  to  overlay  the  exact  solution  in  the  (a)  subplot.  The  (b)  subplot,  therefore, 
displays  the  absolute  error,  i.e.  the  absolute  value  of  the  difference  between  the  numerical  and  exact  solutions.  It  can 
clearly  be  seen  in  this  figure  that  Crank-Nicolson  has  about  an  order  and  a  half  higher  error  than  the  three  ESDIRK 
schemes,  which  themselves  have  quite  similar  error  curves.  This  is  the  expected  result  taking  into  account  the  trends 
seen  in  the  previous  figure. 


X-Coordinate 


X-Coordinate 


(a) 


(b) 


Figure  17:  (a)  Density  and  (b)  absolute  error  in  density  at  the  vortex  center  line  for  the  time  discretization  given  in  the 
legend  with  eighty-one  points  across  the  vortex  at  a  CFL  number  of  about  one  after  convecting  four  periods. 


Figure  (18)  shows  the  density  profile  at  the  center-line  for  the  isentropic  vortex  at  a  CFL  number  of  about  eight  on 
a  grid  with  one-hundred  points  in  the  x-  and  y-dircctions  after  having  convected  for  two  and  four  domain  lengths.  An 
additional  curve  has  been  added  to  these  plots,  designated  ’’CN2-MAX”  corresponding  to  the  profile  at  the  maximum 
density  perturbation  for  the  Crank-Nicolson  time  scheme.  This  additional  curve  was  added  because  the  vortex  has 
meandered  off  the  center-line  by  six  grid  spacings  after  two  domain  lengths  and  by  ten  grid  spacings  after  four  domain 
lengths  of  travel.  As  can  be  seen,  differences  among  all  schemes  besides  ESDIRK4  and  ESDIRK5  are  apparent  at  this 
larger  CFL  number.  Additionally,  it  can  be  observed  that,  although  Crank-Nicolson  does  an  excellent  job  at  preserving 
the  magnitude  of  the  perturbation  (low  dissipation),  the  fact  that,  because  of  dispersion  error,  the  vortex  travels  along 
an  incorrect  path  (it  convects  both  faster  and  at  a  slightly  downward  angle  compared  to  how  it  should  convect),  actually 
leads  to  this  scheme’s  having  higher  error  than  the  more  dissipative,  but  less  dispersive  ESDIRK3  scheme. 
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Figure  18:  Density  at  the  vortex  center  line  for  the  time  discretization  given  in  the  legend  with  eleven  points  across 
the  vortex  at  a  CFL  number  of  about  eight  after  convecting  (a)  two  periods  and  (b)  four  periods. 


Finally,  Figure  (19)  is  analogous  to  Figure  (16),  plotting  temporal  scheme  error  at  a  constant  CFL  number  of  about 
eight  versus  temporal  and  spatial  resolution  for  the  three  ESDIRK  schemes.  As  can  be  seen  from  the  (a)  subfigure, 
ESDIRK3  converges  as  third  order  for  all  resolutions,  having  much  higher  error  at  fine  resolutions.  As  such,  its 
error  dwarfs  the  error  of  the  other  two  ESDIRK  schemes;  thus,  it  was  not  included  in  subfigure  (b).  The  second  plot 
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demonstrates  that  for  ESDIRK5,  error  is  roughly  independent  of  resolution,  while  relative  error  increases  at  better 
resolutions  when  the  ESDIRK4  scheme  is  used.  This  is  the  expected  result  and  points  toward  a  preference  of  having 
the  order  of  accuracy  of  the  selected  time  scheme  be  the  same  as  that  of  the  selected  spatial  scheme  to  guarantee  that 
overall  error  always  converges  at  a  constant  slope. 


dx  (log  scale) 


dx  (log  scale) 


(a) 


(b) 


Figure  19:  Comparison  of  the  density  error  for  ESDIRK4  and  ESDIRK5  time  discretizations  and  a  CFL  number  of 
about  eight  for  (a)  absolute  error  and  (b)  error  scaled  by  the  ESDIRK5  error  at  a  CFL  number  close  to  one. 


V.  Conclusions  and  Future  Work 

The  von  Neumann  analysis  presented  herein  highlights  advantages  and  disadvantages  of  the  various  time  integra¬ 
tors.  As  expected,  Crank-Nicolson,  with  its  low  error  constant  (the  lowest  of  all  second-order  ,4-stable  linear  multistep 
methods),3  performs  well  in  terms  of  dissipation  error,  but  is  more  dispersive  because  of  its  second-order  accuracy. 
Moreover,  its  lack  of  L-stability,  as  is  frequently  cited,4  remains  a  concern.  As  such,  the  six-stage  (utilizing  five  im¬ 
plicit  stages  that  require  solution)  ESDIRK4  appears  to  strike  the  best  compromise  between  diffusion  and  dispersion 
errors,  while  its  stage-count  makes  it  a  candidate,  in  terms  of  efficiency,  for  the  preferred  scheme.  The  compromise 
nature  of  this  choice  becomes  especially  apparent  at  CFL  numbers  greater  than  one. 

The  computational  results  tend  to  support  the  choice  of  ESDIRK4  as  the  preferred  temporal  integrator  for  fifth- 
order  spatial  schemes.  It  performs  just  as  well  as  ESDIRK5  for  typical  time-step  sizes  and  mesh  resolutions.  That 
being  said,  having  the  order  of  accuracy  of  the  spatial  and  temporal  schemes  being  the  same  does  have  the  advantage 
of  guaranteeing  that  the  ratio  of  spatial  to  temporal  error  remains  constant  no  matter  the  grid  resolution  at  a  given, 
constant  CFL  number. 

This  work  clearly  demonstrates  that  the  common  practice  of  using  a  third-order  time  integrator  (like  Runge-Kutta- 
Wray)  coupled  with  a  fifth-order  spatial  discretization  is  inadequate.  As  the  grid  is  refined,  the  time  step  must  be 
refined  to  a  greater  degree  to  experience  an  overall  fifth  order  accurate  decrease  in  the  overall  error.  Specifically,  if 
spatial  resolution  is  doubled,  spatial  error  will  theoretically  decrease  by  a  factor  of  thirty-two  while  temporal  error 
would  only  decrease  by  a  factor  of  eight  if  the  CFL  number  were  kept  constant.  To  get  the  same  factor  of  thirty-two 
decrease  temporal  resolution  would  need  to  increase  by  a  factor  of  25/3  «  3.175. 

Future  work  should  be  divided  into  three  areas.  First,  /,-stable  Runge-Kutta  time  integrators  whose  growth  factors 
decrease  more  slowly  away  from  the  origin  in  the  complex  left  half  of  the  eigenvalue  plane  should  be  developed 
for  use  with  high-Reynolds  number  flows.  Schemes  with  these  properties  (if  they  exist)  should  prove  to  have  lower 
dissipation  error  while  retaining  the  dispersion  error  characteristics  of  high-order  temporal  schemes  already  observed. 
Additionally,  the  present  work  highlights  the  need  for  lower  error  spatial  schemes  with  the  same  formal  order  of 
accuracy.  Such  spatial  discretizations  would  take  maximum  advantage  of  the  high-order  time  discretizations  both 
on  coarser  meshes  and  when  using  CFL  numbers  around  one.  Candidates  for  lower-error  spatial  schemes  include 
compact-difference  schemes 1 6  for  their  enhanced  spatial  resolution  and  filtering  schemes 1 7- 1 8  because  of  their  potential 
for  highly  scale-discriminant  damping.  Finally,  a  third  area  for  future  work  is  to  add  the  preconditioning  techniques 
described  above8  to  the  current  methods.  By  adding  preconditioning,  more  efficient  and  accurate  solutions  should  be 
produced. 
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A.  Butcher  Tableaux 


This  appendix  presents  the  Butcher  Tableaux  for  all  the  of  the  Runge-Kutta  schemes  presented  herein.  The  tableaux 
for  the  non-ESDIRK  schemes  can  be  referenced  in  Butcher19  while  the  tableaux  for  the  ESDIRK  schemes  are  credited 
to  these  sources.1,2 
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